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We investigate stochastic models of particles entering a channel with a random time distribution. 
When the number of particles present in the channel exceeds a critical value N, a blockage occurs and 
the particle flux is definitively interrupted. By introducing an integral representation of the n particle 
survival probabilities, we obtain exact expressions for the survival probability, the distribution of 
the number of particles that pass before failure, the instantaneous flux of exiting particle and their 
time correlation. We generalize previous results for N = 2 to an arbitrary distribution of entry 
times and obtain new, exact solutions for N = 3 for a Poisson distribution and partial results for 
N > 4. 
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I. INTRODUCTION 

A stream of particles flowing through a channel may 
be slowed or blocked if the number of particles present 
exceeds the carrying capacity of the channel. This phe¬ 
nomenon is widespread and spans a range of lengthscales. 
Typical examples include vehicular and pedestrian traf¬ 
fic flow, filtration of particulate suspensions and the flow 
of macromolecules through micro- or nano- channels. A 
specific example of the first category is a bridge that col¬ 
lapses if combined weight of the vehicular traffic exceeds 
a threshold. In filtration, experimental data of the frac¬ 
tion of grains retained by a filter mesh can be explained 
by assuming that clogging may occur when two or more 
grains are simultaneously present in the same vicinity of a 
mesh hole, even though isolated grains are small enough 
to pass through the holes [1]. A biological example is 
provided by the bidirectional traffic in narrow channels 
between the nuclear membrane and the cytoplasm[2]. 

The totally asymmetric simple exclusion effect process 
(TASEP) provides a theoretical approach to these phe¬ 
nomena. The TASEP is a lattice model with a stochastic 
dynamics where particles hop randomly from site to site 
in one direction with the condition that two particles 
cannot occupy the same site at the same time [3, 4]. At 
the two extremities of the finite lattice, particles are in¬ 
serted and removed with two different rates. The model 
and its extensions provide quantitative descriptions of the 
circulation of cars and pedestrians[5-10]. The so-called 
bridge models consider two TASEP processes with oppo¬ 
sitely directed flows, but allow exchange of particles on 
the bridge[ll-15]. At the microscopic level active motor 
protein transport on the cytoskeleton has been modeled 
by a TASEP [16, 17]. 

Recently, some of the present authors [18, 19] intro¬ 
duced a class of continuous time and space stochastic 
models that are complementary to the TASEP approach. 
In these models particles enter a passage at random times 


according to a given distribution. In the simplest concur¬ 
rent model particles move in the same direction and an 
isolated particle exits after a transit time r but if N = 2 
particles are simultaneously present, blockage occurs. If 
the particle entries follow a homogeneous Poisson process 
all properties of interest, including the survival probabil¬ 
ity, mean survival time and the flux and distribution of 
exiting particles can be obtained analytically. The model 
has a connection to queuing theory in that it is a gen¬ 
eralization of an M/D/1 queue, i.e. one where arrivals 
occur according to a Poisson process, service times are 
deterministic and with one server. This queue has many 
other applications including, for example, trunked mobile 
radio systems and airline hubs [20-22]. 

Opposing streams, where blockage is triggered by the 
simultaneous presence of two particles moving in different 
directions can be treated within the same framework [18]. 
Inhomogeneous distributions of entering particles can be 
treated analytically [23]. It is also possible to obtain 
exact solutions for when the blockage is of finite duration, 
rather than permanent [24]. In this case, for a constant 
flux of incoming particles the system reaches a steady- 
state with a finite flux of exiting particles that depends 
on the blockage time it,. 

The purpose of this article is to explore the properties 
of the concurrent flow models for any distribution of en¬ 
try times and when the threshold for blocking is IV > 2. 
In addition to the applications described above, this gen¬ 
eralized model may also be relevant for internet attacks, 
in particular denial of service attacks (DoS) and a dis¬ 
tributed denial of service attacks (DDoS) where criminals 
attempt to flood a network to prevent its operation [25— 
27]. 

Unfortunately, the method used to solve the models for 
N = 2 [18, 24] applies only to a Poisson distribution and 
cannot be used even in this case for N > 2. In section 
II, we develop a new approach providing formal exact 
expressions of the key quantities describing the kinetics 
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FIG. 1. Concurrent flow model N = 3. Particles enter the 
left hand side of a channel of length L randomly. Top: two 
particles cross and exit the channel in a time r. Bottom: If a 
third particle enters while the two previous particles are still 
in the channel, a blockage occurs instantaneously. 


of the model. In section III, as a first application, we 
recover the results of the model N = 2 that were first 
obtained by using a differential equation approach[18]. 
In section IV, we present a complete solution when the 
entry time distribution is Poisson for N = 3. In section 
V we consider the case of general N. In section VI we 
investigate the time correlation for N = 2 and N = 3, 
and we further explore the model by studying the corre¬ 
lations between the arrival times of the particles. We also 
explore the connection with the equilibrium properties of 
the hard rod fluid. 


II. CONCURRENT FLOW MODEL 
A. Definition 

We assume that at t = 0 the channel of length L is 
empty. The first particle enters at a time t 0 that is dis¬ 
tributed according to a probability density function ip(s). 
The entry of subsequent particles is characterized by the 
inter-particle time ti,i > 0 between the entry of parti¬ 
cle i and i + 1. We assume that the t, ; are distributed 
according to ip(s) and uncorrelated. The total elapsed 
time is then t = to + Y'™ti +t' when n particles have 
entered and t' is the time elapsed after the entry of the 
last particle. 

If unimpeded by the presence of another particle, a 
particle exits after a transit time r > 0. Blockage occurs 
when N particles are present in the channel at the same 
time, which occurs if ti + f,;+i + • • • + tj+jv -2 < t (see 
Fig.l for the case N = 3). The model is non-Markovian 
as the state of the system at time t depends not only on 
the actual state but also on the history of the system. 
The probability that no particle enters in the interval 
[0,f] is 1 — ip c (t) with ip c {t) the cumulative distribution 
^ c(t ) = /o i/j(s)ds. 

The simplest case is a homogeneous Poisson process 
where the probability density function of particle times 
is ip{t) = Ae~ At where A is the rate (sometimes called the 
intensity). 


B. Quantities of interest 

The key quantities describing the process are the prob¬ 
ability that the channel is active at time f, namely the 
survival probability, p s {t ), the average blocking time ( t) 
(where the bracket indicates an average over realizations 
of the process), the number of particles that have ex¬ 
ited the channel at time t, ( m(t )), and the instantaneous 
particle flux j(t). 

The survival probability can be expressed as the sum 
over all n-particle survival probabilities q{n,t) 1 i.e. the 
joint probability of surviving up to t and that n particles 
have entered the passage during this time, 

OO 

Ps{t) = ^2q(n,t) (1) 

n—0 

For general N and n > N — 1, q{n,t) can be expressed 
as: 


q{n,t) = / 
Jo 


]Q dti^(t. 

i—0 

n-AT+1 / N-2 

n 0 ( ^+ m _ r 

j =1 \m —0 


dt’{ 1 - ip c (t')) 

/ n—1 

* t-5> 


1=0 


t’ 
( 2 ) 


where 0(x) the Heaviside step function. The first n inte¬ 
grals correspond to the arrival of n particles in the chan¬ 
nel, with time intervals tj, the integral over t' imposes 
that no particle enters after particle n. The Heaviside 
functions account for the constraint that no consecutive 
sequence of N particles can be simultaneously in the 
channel, i.e. in a time interval smaller than r and the 
S function imposes that the observation time t is equal 
to the sum of the time intervals ti plus to and t'. 

For 0 < n < N — 1 there is no constraint on the par¬ 
ticle time interval so the probability q(n,t) is expressed 
as the joint probability of n independent and identically 
distributed events 


q{n,t) = / 
Jo 


‘n—1 


dti^U) 

2—0 

n—1 \ 

5 (t ^ ti -1' j 


dt'( 1 - ipc (*')) 


( 3 ) 


i =o 


and 


g(0,f) = 1 - ip c {t) 


( 4 ) 


Once the q(n,t), and hence p s (t ), are known we can 
obtain several useful quantities. The probability density 
function of the blocking time, f(t) is simply related to 
Ps(t) 


m 


dp s {t) 
dt 


( 5 ) 
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Defining the Laplace transform as f(u) = / 0 °° dte ut f(t), 
one infers 

f(u) = 1 - up s (u) ( 6 ) 

The mean blocking time is given by 

POO 

(t) = / dttf{t)=p s { 0) (7) 

J o 

The instantaneous flux of particles exiting the channel 
can be obtained by noting that if a particle exits the 
channel at time t, at most N — 1 particles can enter the 
channel between t and t — r if no blockage is to occur. 
Since blockage is irreversible the flux tends to 0 when the 
time increases, j( oo) = 0 for all value of N. The total 
flux is given by the sum, 

OO 

j{t) = ^2j(n,t) ( 8 ) 

n— 1 

where j(n, t) is the partial flux where a particle exits the 
channel at time t such that the channel is still open and 
n particles have already entered, for n > N 

'n— 1 

.2=0 

n-iV+1 /N-2 

II M 55 tj+rri — T 

j =1 \m—0 

N—2 k 

[s(t r -r) + y w + 55 tn ~ w - r ))]> n > N 

k= 1 w—1 

( 9 ) 

The condition that a particle exits at time t is ex¬ 
pressed in terms of <5 functions. More specifically the 
exiting particle can be the last particle to enter, corre¬ 
sponding to the term S(t' — r), or one of the other TV — 1 
previously entering particles, corresponding to the sum 
over 5 functions. For n < N blocking is not possible, 
so Eq.(9) is replaced by one without the Heaviside func¬ 
tions. 

Finally, the number of particles that have exited at 
time t can be obtained by integrating over the particle 
flux 

(m(t))= f dt'j{f) (10) 

Jo 

We can also obtain the distribution of particles exit¬ 
ing the channel. Let h(m, t) denote the probability that 
blockage occurs in the interval ( 0 , t) and that m particles 
have exited during this time. Its time evolution is given 

by 


s (t - -1' 


i =0 


J'O 


ht)=r 

Jo 


dt'{ 1 - V’c(f')) 


dh(m, t) 
dt 


oo m-\-N—l m / N—2 

dtii>{ti) n 0 t -i+p ~ t 

i—0 j—1 \ p—0 

N— 1 m-\-N— 1 

0(t - Y, t m+p )S(t- ti)m> 1 (11) 

p= 1 i=0 



The upper part of the right hand side corresponds to 
the event where m + N particles have entered at time t, 
and there was no blockage involving the first m + N — 1 
particles. The second Heaviside function corresponds to 
the constraint that the last N particles are blocked in the 
channel, with the N + m th particle entering at time t. 

One can check that 


OO 

( m(t )) = mh(m.,t) 

m =0 



( 12 ) 


We now consider the specific cases N = 2 and N = 3. 


III. N = 2 


Since each Heaviside function in Eq.(2) depends on 
only one variable, the multiple integrals can be always 
calculated. Taking the Laplace transforms of Eq.(2) and 
Eq.(3), one obtains 


q{n,u) = ip{u)( - if> c (u)) 

u 


dte 


(13) 


Using Eq.(l) and ^ c (u) = we obtain the Laplace 
transform of the survival probability. 


Ps(u) = yq(n,u) 


n—0 


1 — V’(w) 


i + 


%i{u) 


1 — f T e~ ut ip{t)dt ^ 
Therefore, the mean time of blockage is 

1 


(t) = t 


1 + f 0 T ^) dt 


(14) 


(15) 


where t = ip'(0) = f£° dtti/)(t) is mean inter parti¬ 
cle time. To interpret Eq.(15) we note that ^> c (r) = 
f 0 ip(t)dt gives the probability that two consecutive par¬ 
ticles are separated by a time smaller than r. 

For a Gamma distribution, ^(t) = \ a t a ^ l e~ xt /T{a) 
where a is a shape parameter, the mean time of blocking 
is equal to 


<*> 


a 

A 



r(q) \ 

F(a) - 7 (t, a)J 


(16) 


where r(a) and 7 ( 0 :, x) are the Gamma and incomplete 
Gamma functions, respectively. When Ar < 1, one ob¬ 
tains 


<*> 


1 a! 

A (Ar) a 


(17) 


Figure 2 shows < t > versus Ar for a = 2,3,4. One ob¬ 
serves an excellent agreement between simulation data 
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FIG. 2. Mean time of blocking < t > as a function of the in¬ 
tensity, A, for a Gamma distribution for a = 2, 3,4 (from bot¬ 
tom to top), from numerical simulation (circles) and Eq.(15) 
(full lines). Dotted lines correspond to the asymptotic behav¬ 
ior, Eq.(16). 


FIG. 3. Mean flux j(t) as a function of time for a Gamma 
distribution with a = 3 and for different values of A (r = 1). 
The solid lines show the exact expression, Eq.(21) and the 
circles show simulation results. 

In Laplace space, the summation over n can be per¬ 
formed and j ( u ) is given by 


(circles) and the exact formula, Eq.(16). As expected, 
the mean time < t > diverges when A goes to zero. The 
asymptotic behavior, Eq.(17) provides a good approxi¬ 
mation of simulation data when At < 1. 

By taking a = 1 in the Gamma distribution, which 
corresponds to a homogeneous Poisson process, the mean 
time of blockage is given by 


(t) = 


2 — e 


— A r 


A(1 - e~ Xr ) 


(18) 


a result previously obtained by using a master equation 
for the time evolution of the q(n,t) [18]. 

The mean flux j(t) can be obtained by using Eqs.(8,9) 


°° /»00 /* C 

„_1 Jo Jo 


n—l 
n— 1 




i=0 


n e - t ) 

t=i 


n—l 


$ t - E ~ M - 


i=0 


(19) 


The multiple integral can be factorized and the flux is 
given by : 


j(t) = 


n—l * 


dt 0 ip(t 0 ) / dt'(l - 


d ^ - E 



( 20 ) 


j(u) = 


(1 - i>c(r))e UT '0(u) 






With a Poisson distribution = Xe At , we have 

Ae -(«+A)r 


i(w) = 


, + A(1 - e-(“+ A ) T ) 


( 21 ) 


( 22 ) 


By taking the inverse Laplace transform, the mean flux 
j{t) can be expressed as a series 


j(t) = Ae- A 'E 

n—l 


A(A it - (n + l)r))"6»(A(f - [n + l)r)) 
n\ 


(23) 

as obtained previously by using a master equation ap¬ 
proach [19]. 

No particle exits the channel between 0 and r; indeed, 
the flux is obviously equal to 0 in this interval and rises 
instantaneously to a maximum, j max = Ae _Al " which it¬ 
self is maximum when A = and then decreases to 0. 

For a Gamma distribution with an integer value of a, 
the Laplace transform of the flux can be obtained explic¬ 
itly, but increasing a it rapidly leads to lengthy expres¬ 
sions. 

Figure 3 displays the time evolution of the mean flux 
j{t) for different values of A and a = 3. In all cases, 
the flux becomes nonzero for t > r, corresponding to the 
exit of a first particle. For At > 1, j(t) displays a strong 
maximum at a a time t m slightly larger than t and decays 
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to 0. For At = 0.5, the maximum of the flux is shifted 
to a time t m — 3r and the typical decay time is around 
lOOr. For At = 0.25, j(t') increases up to a quasi-plateau 
and the typical decay time is larger than IOOOt, which 
corresponds to a physical situation where a large number 
of particles exit the channel before the definitive clogging. 
Note that for a given value of A the flux is much larger 
than for a Poisson distribution. However, it approaches 
zero for sufficiently long times with a characteristic time 
equal to the mean blocking time. 

We also consider the probability, h(m, t), that blockage 
occurs in the interval (0, t) and that during this time m 
particles exit the channel. The time evolution of this 
function is given by 


dh(m, t) 
dt 


«oo m+1 m 

= / n dtiip(ti) e (tj - t) 

"'° i =0 j=l 

ra +1 

6{t - - ^2 *0 (24) 

i =0 


Two particles have to be in the channel for the system 
to block, so the interval between 2 consecutive particles 
has to be less than t (the 6 function). The previously 
entering particles exited the channel without blockage. 

Taking the Laplace transform we obtain for m > 0 


h(m, u) 


$(u) 

u 




—ut 


dt 1 



dtip(t)e ut 

(25) 


The probability that the channel is blocked can be ex¬ 
pressed as the sum over partial probabilities h(n 7 t ), 
namely h(t) = oh{m,t). By using Eq.(25), one 
infers lim h(t) = lim uh(u) = 1, as because block- 

t—¥ OO U— >-0 

age is certain to occur, a result valid for any distri¬ 
bution ip(t). Finally, we note the following sum rule, 
Etj>ofe( n ,^) + h(n 7 t)) = 1 - all configurations of the 
process are either blocked or unblocked. 

For the Poisson process, an explicit expression can be 
obtained 


h.(m, w) = - 


vm+2 


(A + u ) 


m+2 


l _ 6 -( a +“)t 


0 —(A+ii)mr 


(26) 


Performing the Laplace inversion we obtain h(m 7 1 ) as 
obtained previously [19]. As expected, h(m,t) is equal 
to zero for t < mr corresponding to the minimum time 
necessary for m particles to exit the channel. For the 
Gamma distribution with a = 2 we obtain 



FIG. 4. The probability h(m, t ) as a function of time for a 
Gamma distribution with a = 2, A = 2 and m = 0,1, 2 (from 
top to bottom). The full curves show the exact expression, 
Eq. (27) and circles show simulation results. 


for a = 2 and A = 2. As expected, h(m , t) = 0 for t < mr 
, which can be explained by the fact that the minimum 
time for having a configuration where m particles exit 
the channel must be at least larger than mr. Similarly, 
the transient time associated with him, t ) increases with 
to, and corresponds to rare events when to increases. 


IV. N = 3 


For the first three partial probabilities, there is no con¬ 
straint and one easily obtains that <?(0,f) = l — ip c (t)i 
and for * = 1,2 the probabilities are given in terms of 

the Laplace transforms q(i,u) = u) 1 . For 

a Poisson process, one recovers that g(0,t) = e ~ xt , 

q(l,t) = A te~ M and q(2,t) = For n > 2, 

Eq.(2) becomes 


/*oo /*oo n—1 

q(n,t)= dt'{\ - i!> c (t')) / TTdfiV^i) 

Jo Jo i=0 


n —2 


n— 1 


x n + t J+ i - t - x! ^ ~ o ( 2s ) 

3 =1 *=0 


h(m,u) = 


/g— (»+A)rm^ 2 (m+ 2 ) 
I ^ 


u(u + A) 2 ( m + 2 ) 

(1 + {u + A)t)” 1 (1 - e~ T ^ u+x \l + (u + A)t))) 

(27) 

For the Gamma distribution, we plot in Fig. 4 the time 
evolution of h(m 7 1 ) as a function of time with to = 0 ,1, 2 


The constraint, imposed by the 6 function, requires that 
the sum of two consecutive time intervals be less than t. 
Taking the Laplace transform of Eq.(28), one obtains 

{(„,») = MeriMl 1*° dt*p(t)e~ ut r(n - l,t,u) 

U Jo 

( 29 ) 
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where the auxiliary function r(n — l,t,u) is given by 

„oo n-2 n —2 

r{n -l,t,u)= dti%lj(ti)e~ uti 9(tj + t j+1 - r) 

•'0 »=1 j —1 


(30) 

where t n — 1 = t. A recurrence relation can be written for 
r(n, t, u ) 

nOO 

r(n,t,u) = / dt'ip{t')e~ ut r(n — (31) 

J max(r—i,0) 


B{z,u) are determined by Eq.(35) adapted to a Poisson 
process. 

For n = 0,1,2 the partial probabilities q(n,t) corre¬ 
spond to those of a Poisson process. For n > 2, by using 
the generating function G r (z,u), the Laplace transform 
of q(n, t ) is given by 


q{n,u) 


(X + u) 2 


dt\e~ {x+u)t 


d n 2 G r (z,t,u) 
dz n ~ 2 


z=0 

(38) 


After some calculation, one obtains 


with r(l, t, u) = 1. 

Let us introduce the generating function G r (z, t , u) de¬ 
fined as 


G r (z,t,u) = ^ 1 r(n,t,u) (32) 

n=l 

Multiplying Eq.(31) by z" -1 and summing over n, one 
obtains that 

/»oo 

G>(, 2 , t, u) = 1 + ^ / difip(t')e~ ut G r (z, t', u) 

J max(r—1,0) 

(33) 

For t > t G r (z,t,u) is constant, i.e. G r (z,t,u) = 
G r (z,r,u). For t < t, it is convenient to express the 
time evolution of G r (z,t,u) as follows: taking the first 
two partial derivatives of G(z, t, u) with respect to t, one 
obtains the ordinary differential equation 

d 2 G r (z,t,u) ( ip(T — t) \ dG r (z,t,u ) 

dt 2 y ip(T — t) + U J Qt 

— z 2 ip(r — t)ip(t)e~ UT G r (z, t) (34) 

By using Eq.(33), the differential is supplemented by two 
boundary conditions 


<z(3, t) =0{t-T) A 3 e xt [kr(t - t) 2 + ±(t - t) 3 ] 

<7(4, t) =A 4 e- At \e{t-T) { ^-9{t-2T)^=^ 

(39) 

By using Eqs.(29) and (32), the Laplace transform of 
the survival probability p s (t) is 


p s (u) =q{ 0, u) + q{ 1, u)+ 


dtip(t)e wt G r (l,t,u) (40) 


By inserting the solution of Eq. (34), the Laplace trans¬ 
form of the survival probability is given by 


Pa(u) = 


A 


(A + u) 2 
+ B(l,u) 


1 + — + A{ 1, u) 
1 + — (1 — e~ SlT ) 

Sl 


1 + —(1 - e 
S 2 


2T ) 

(41) 


where 


A(l,u) = 
B(l,u) = 


Ae S2T (s 2 - A)(si + s 2 ) 
A 

Ae SlT (si - A)(si + s 2 ) 


(42) 


G r (z, 0 ,u) 

dG r (z,t,u) 
dt 


= 1 + zG r (z, r, u) dt'ip(t')e ut 
= zip(0)G r (z, 0, u ) 

(35) 


Eq.(34) cannot be solved analytically in general but 
for a Poisson distribution it becomes 

d 2 G r (z, t, u) dG r (z,t,u) 

—ae— = (A + “ ) — at — 

- (z\) 2 e~( u+x)T G r (z,t,u) (36) 

with the boundary condition given by Eq.(35) with 
ip(t) = e~ xt . 

The solutions of the characteristic equation of Eq.(36) 
are 

, , (A + u)± \J (A + u) 2 - 4 (zX) 2 e~A+n)r 

Sl,2(Z,U) — - - - 

(37) 

and the generating function is given by G r (z,t, u) = 
A(z,u)e s A z < u ) t + B(z,u)e s A z ' u ') t where A(z,u) and 


with 

A = e (si+S2T) sis 2 (si - s 2 ) + A (s 2 2 e S2T - s?e SlT ) (43) 


From the generating function, one can also obtain 
global quantities, like the mean blocking time (t) = p s (0 ) 
Let g = y/| 1 — 4e _AT | and v = then, after some 
calculation, one obtains for Ar > 2 ln(2) 


A (t) 


2e u sinh(^) + ge Xr 

—g — 2 sinh(^)e -1 ' + e v (sinh(zv) + g cosh(^)) 


and for Ar < 21n(2) 


+ 1 
(44) 


A(t> =- , 2e" sinjO + y Ar - +1 

—g — 2 sm(^)e _1/ + e v (sm(i/j + g cos(is)) 

(45) 

Fig. 5 shows the mean blocking time (t) of the mod¬ 
els with N = 2,3,4, 5, 6, 7 for a Poisson distribution ob¬ 
tained by simulation and for N = 2,3 by using the an¬ 
alytic expressions. We observe a perfect agreement be¬ 
tween simulation data and exact expressions for N = 2 
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FIG. 5. Mean time of blocking as a function of the intensity At 
for IV = 5,4, 3, 2, top-to-bottom, from numerical simulation 
(circles) and Eq.(18) N = 2 and Eqs.(44,45) N = 3 (full 
curves) for a Poisson distribution. The inset compares the 
asymptotic formula, Eq.(58), with simulation results. 


Eq.(18) and N = 3 Eq.(44,45). More generally, one 
observes a divergence of the mean blocking time as At 
goes to 0 and indeed performing a first-order expansion 
of Eq.(45) in At gives 


FIG. 6. Mean flux j(t) as a function of time for a Poisson dis¬ 
tribution for A=1,2(t = 1). The solid lines show the exact 
expression, (inverse Laplace transform of Eq.(49)) accurately 
matching simulation results (wavy lines). 


be expressed as 


«> = jlp <«> 

The mean flux j(t) can be also obtained by using 
Eq.(9) and the auxiliary functions r(n,t,u) and it comes 
for the Laplace transform j(n,u) (for n > 1) 


iO) 


Ae -Cu+A)r 

A + u 


- 4 ( 1 , u) 




+ 13(1,u) 



X_ 

si 


(49) 


where -4(1,it) and 73(1, it) are given by Eq.(42). 


j(n,u)=e “ T t/((u) ^(1 - V’c(t)) J dte ut i!){t)r{n - 1, t, 
+ J dtij}(t){l - tp c (T - t))r(n - 

(47) 

By summing over n (accounting for the boundary 
terms j(l,it) and j(2,u), the Laplace transform j(u) is 
expressed as 


Because the right-hand-side of Eq.(49) can be factor- 
' ized by e~ UT , it implies that j[t) = 0 for t < r, which 
corresponds to the minimum time for a particle to exit 
the channel. 

The mean flux j(t) is plotted as as function of time 
for A = 1, 2 with a Poisson distribution for A = 1, 2 (Fig. 
6 ). A discontinuity appears at t = t where the flux is 
maximum j(r) = A. At t = r, the flux is given by 

j( T ) = A(! + Xr)e~ XT (50) 


j(u) =e “ r V)(u)(l - V’c(t)) / dte ut i/)(t)G r (l,t,u) 

Jo 

-f e~ UT ^){u) f dt^(t)(l — ^ c (t — t))G r (l,t, u) 

Jo 

+ jO,u) (48) 

By using Eq.(35) and the expression of the generating 
function G r ( 1, t, u), the Laplace transform of the flux can 


which corresponds to events where a particle exits be¬ 
tween t and t + dt such that 0 or 1 particle is still in 
the channel. The flux decay exhibits a visible cusp at 
t = 2 t which corresponds to the non analytical structure 
of the solution. At long times, the flux decays to 0, with 
a typical time which becomes larger when A decreases. 

The joint probability /i(m, t ) can also be obtained with 
the function r(n,t,u). For m > 1 its time evolution is 













V. TV > 4 



FIG. 7. Probability distributions h(m , t) versus time t for a 
Poisson distribution, for m = 0,1, 2 (from top to bottom) and 
A = 1. The solid lines correspond to the model with TV = 2 
and dashed lines correspond to the model with TV = 3. 


given by 

77 / ,n r oo ™+2 m 

-ijr = / n n 0 + *j+i ~ r ) 

70 j=0 J=1 

ra +2 

0(t - tm + 1 - ^ ti) (51) 

2=0 


Taking the Laplace transform gives 

7/ \ /*oo 

h(m,u ) / TT dtiip(ti)e~ uti 

u J o i= 1 

m 

JJ 0(tj + t j+ 1 - t)8(t - t m +1 - i m+2 ) (52) 

i=i 

that can be expressed using the function r(n,t,u ) as 

h(m,u) ( dti/j(t)e~ ut 

u Jo 

r T 

/ dlJip(l/)e ut r{m + l,t',u) (53) 

Jo 

Figure 7 shows the time evolution of the probability 
distributions h(m,t) for a Poisson distribution with A = 
1. The probability that zero or one particle (m = 0,1) 
particle exits is smaller for N = 2 than for = 3. For 
m > 2 the order reverses (e.g. for m = 5, case shown). 
This is because, for a given value of At, more particles 
exit before blockage as N increases. 


We have seen that for TV = 3 the product of Heavi¬ 
side functions in Eq. (2) leads to a simple recurrence rela¬ 
tion Eq.(31). For N > 4 the task is much more difficult 
because one needs to introduce auxiliary functions that 
depend on N — 2 time variables. These functions are re¬ 
lated by an integral equation that cannot be converted 
to an ordinary differential equation. We therefore pro¬ 
pose an approximate treatment of the dynamics. For the 
model where the blockage occurs when N particles enter 
the channel between t — r and t, the first TV — 2 partial 
probabilities q(i, t) obey differential equations identical 
to those of a Poisson process 

^=—A,(0 ,*) (54) 


and 

^ ^ = —A q(n, t) + \q{n —1,2), 1 < n < N — 1 (55) 

dt 

For n > N — 1, the non-Markovian constraint applies, 
but for n = TV, the time evolution is simply given by 


= ~X„(N,t) + X N f bdV-,«V-l-M-r) 


8 =0 


(56) 

The gain term reflects the fact that blockage only oc¬ 
curs with TV particles, TV — 1 terms correspond to the 
cases where there may be from 0 to TV — 1 particles in 
the channel. 

For n > TV, the dynamics of q(n, t) for n > TV can be 
approximated as follows 


dq(n, t) _ ^ _j_ \q^ n — 1 , f _ T )e Xt 

dt 

N—2 T 

+ AV / dt 1 K s (ti)e~ XT q s (n - 1 — s,t - r - ti) 

3 = 1 

(57) 


where we have introduced a kernel K s (t). We then con¬ 
sider two physical situations. In the first, the last s 
particles are assumed to have entered the channel in 
an infinitesimal time interval and the kernel is given by 
K s (t) = This choice overestimates the survival 

probability. TV — 2 particles can be in the channel (so can 
enter between t — r and t ) when a new particle enters. 
The other particles enter between time 0 and t — r. This 
fails to take into account some blocking. In the second 

case we take K s (t ) = e ~ Xt which is proportional 

to the probability that s — 1 particles enter in (0,2). This 
choice underestimates the survival probability. When a 
particle enters at time t there may be a maximum of TV—2 
particles in the channel to avoid blocking. If a particle 
arrives at a time 2 i between t — r and 2 , there may be a 
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maximum of TV — 3 particles between t and t — t\ and no 
particle between t — t\ and t — t\ — r. 

Taking the Laplace transform of Eq.(57), we calculate 
two different generating functions corresponding to the 
two kernels, and the corresponding mean survival times. 
These bracket the exact value and for Ar <C 1 the two 
solutions approach the same limit: 


(t) 


(A r) N 


(58) 


To obtain exact results for TV > 4 is a challenging 
problem. We therefore finish this section by presenting 
some numerical results that illustrate the general trends. 
The inset of Fig. 5 compares the asymptotic behavior 
for mean blocking time, Eq.(58) with simulation results 
for TV = 2 to TV = 5. We observe that the scaling law 
provides provides a good description of the process for 
At < 0.5. 

In Fig. 8 we present numerical results for the mean flux 
of exiting particles as a function of time. This quantity 
acquires a non-zero, maximum, value at t = t given by 


N—2 


j 0 ) = a 

2=1 


M,-- 


(59) 


This expression corresponds to events where a particle 
exits between t and t + dt such that 0,1.., TV — 2 particles 
are still in the channel. For t > t, we observe a drastic 
increase of the characteristic decay time as TV increases 
(see the lower figure of Fig. 8). For TV = 2, j(t) is very 
small for t > 3r, while for TV = 7, the flux is almost 
constant during two decades. 


VI. CORRELATIONS 


We now consider the time correlation function C(t ) 
that represents the density function that any two parti¬ 
cles have a time separation t. C(t ) can be expressed as 
the sum of partial correlation functions c(n, t) that cor¬ 
respond to the probability density that the first and last 
particles of sequence of n + 1 particles are separated by 
t. 

OO 

C(f) = £>,*) (60) 

71=1 

The partial correlation function c(n, t), the joint proba¬ 
bility of having a particle at t = 0 and the nth particle 
at time t , can be written as 


poo iV—2+n iy -4 -til 

c(n, t) = / dtiC (iv_2) (Ti, ...,t N - 2 ) 

1 j=N-l 

( JV+ti-2 \ 

t ^ ^ I 

i=N—l J 

( 61 ) 


iV-2+7i 


= 1 

7i / N—2 

n ^ — t 




FIG. 8. Mean flux j(t) versus time t for TV = 2,3,4, 5,6,7 
(from bottom to top) for a Poisson distribution with A = 1. 
Top: short time behavior. At t = 1, circles correspond to the 
exact values of the mean flux, Eq.(59). Bottom: linear-log 
plot showing the long-time behavior, (r = 1) 


where c^ N 2 \t\, ..., tN- 2 ) is the joint probability of hav¬ 
ing TV — 1 particles such that the first and the second 
particles are separated by a duration of £ 1 , the second 
and the third particles by a duration of and the 

N — 2 and TV — 1 particles by t-N- 2 - We can write this 
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probability as 

c {N - 2) (h,t N _ 2 ) = J dt 0 c (N ~ 2) (t 0 , tjv- 3 ) 

N—2 

x 4>(t N - 2 )0('^2 tj - t) (62) 
3 =1 

This definition of the correlation function considers all 
trajectories, including those that end before a given time 
t. As a result, the correlation function approaches zero 
at long time. It seems more interesting to keep only 
trajectories which have survived until at time t. 

To generate a infinite sequence of particles correspond¬ 
ing to a trajectory of the model, let us consider the fol¬ 
lowing rejection-free algorithm: Accounting for the con¬ 
straints of the model (only less than N particles must 
enter the channel in the duration of time r) without inter¬ 
rupting the traffic, one introduce the discrete stochastic 
equation 


N—2 

t n = max(r - Y t n ~j, 0) + 77 (63) 

3 =1 

where p is a random number generated from the ip dis¬ 
tribution and t n _j,j = 1, N — 2 are the time intervals of 
between the N—2 previously entering particles. 

In order to compute the correlation function associ¬ 
ated with this rejection-free algorithm, we replace ip(U) 
in Eqs. (61,62) with 


N—2 

ip(ti - max(r - Y U-j, 0)). (64) 

i =1 

The partial correlation function c(n, t) can be also ex¬ 
pressed as the average over the event of having a first and 
n + 1 th particles separated by a time duration t 

n 

c(n,t) = (S(t-^2ti)) (65) 

i—1 

The conservation of the probability reads 

dtc(n, t) = 1 ( 66 ) 

By summing over n, the integral correlation function 
C(t) is given at long time by 

dt'C(t’) = (n(t)) (67) 

where ( n(t )) is the mean number of particles along a tra¬ 
jectory for a time duration t. At large t, this quantity 
goes to a constant because we only consider trajectories 
that have survived. By using that C(t) goes to a con¬ 
stant at long time (due to to the decay of the memory 
between particles that entered with a large time differ¬ 
ence) (C(t) —> Coo), we infer that = 1 /t where t is 




the average separation in time between successive parti¬ 
cles. 

We now focus on N = 2 and N = 3 by using the 
rejection-free trajectories for which exact solutions can 
be obtained. 


A. N = 2 


The partial correlation function c(n, t ) is simply given 
as the product of integrals on each independent interval. 
Eq.(63) is very simple t n = r + rj, which means that 
ip(t) is replaced with ip(t — r) in Eq.(64). Therefore, the 
Laplace transform of c(n, t) is given by 


c(n, u) 




c(l,u) n 


( 68 ) 


This results from the fact that successive events are not 
correlated. 

Inserting Eq.( 68 ) in Eq.(60) we obtain 


C(u) 


c(l,u) 

1 — c(l, u) 


(69) 


At long time, C(t) approaches a constant value cor¬ 
responding to a constant mean density. By using the 
factorization property, c(n,u) = c(l,u) n , and the ex¬ 
pansion c(l,u) = c(l, 0 ) + udc(l, u)/du\ u =o + 0(u 2 ) one 
can show that C(oo) = lim^^o uC(u) = 1 /t where 
t = f 0 tc(l,t)dt = —dc(l,u)/du\ u =o is the average in¬ 
terval between particles. That is, the smaller the average 
separation in time between successive particles, the larger 
the steady state value of the time correlation function. 

For a Poisson distribution ip(t) = Xe~ xt we find 


C(u) = y. 


A + u 


(70) 


The inverse Laplace transform gives an explicit expres¬ 
sion 


C(t) = X *Y d(X(t — nr)) 

n =1 


(X(t — nr)) n_ 1 e _A ^ _ ” T ) 
(n - 1 )! 


(71) 


Figure 9(a) shows C(t) for two values of At. As ex¬ 
pected, C(t) is strictly equal to 0 for t < t since no 
particle can be inserted if the delay between two succes¬ 
sive particles is less than r. The maximum of C(t) is 
obtained at t = t where C(t) = A and decreases to 0 at 
large t. Note that a cusp is present at t = 2r, a similar 
behavior observed for the other quantities such as the 
flux and the survival probability. In the long time limit 
C(t = 00 ) = lim^o uC(u ) = A/(l + At). 

It is also interesting to note that correlation function 
Eq.(71) corresponds to the density correlation function 
of the positions of the particle centers in a hard rod fluid 
of density p with A = p/( 1 — p). 
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FIG. 9. Correlation functions for (a.) N = 2 and (b) N = 3 
for A = 1 and 0.5 (lower curves). The solid lines correspond 
to Eq. (71) (N = 2) and Eq. (90) (N = 3) and circles to 
numerical simulations. 


B. N = 3 

For N = 3, the discrete stochastic equation, Eq.(63), 
becomes 


t n = max(r - t„_i,0) + (72) 

where t n denotes the time interval between the n — 1 
and n particles and r] is a random number chosen with 
an exponential probability distribution Xe~ xt . In queu¬ 
ing theory this equation is known as the Lindley-type 
equation [28-30]. 

For the Poisson distribution ip(t), Eqs.(61,62) with 
Eq.(64) gives 


poo n 

c(n,t)= dt 0 c(l,t 0 )6(t 

i= i 


IK 


d£ i Ae _A(t< “ max(r_ti - 1 ’ 0))) 


^/max(r—£i_i,0) 


(73) 


and 

cOO 

c(l,t) = / dtic(l,ti)Ae- A( *- Ma:E(T - tl - 0)) (74) 

J Max(r—t, 0) 

Note that the constraint applies to two consecutive in¬ 
tervals, i.e the arrival time between 3 consecutive parti¬ 
cles is greater than r. Consequently, the partial corre¬ 
lation c(n, u ) is never the product of smaller correlation 
functions, as for the N = 2 model. 

Because the kinetics were obtained exactly in the pre¬ 
vious section only for the Poisson distribution, we restrict 
our analysis to this distribution. 


From Eq.(74), one easily shows that c(l,t) is constant 
for t > t. For t < r, by taking the derivative of Eq.(74), 
one obtains 

dC ^ ^ = A(~c(l, t) + 6{t - t)c{ 1, r - t)) (75) 

whose solution is 

<i, t ) = ~ t) + e ~ Kt ~ T) °- r )) ( 76 ) 


One can easily obtain the average time between two 
consecutively particles in a trajectory. 


t= dtc(l,t)t 

Jo 


(At + l) 2 + 1 
2A(Ar + 1) 


(77) 


As might be expected, when the intensity At is high, 
the probability distribution is uniform within the first 
interval [0, t] and equal to |. Conversely, when At tends 
to 0 the effect of the constraint is negligible, i diverges 
as y, corresponding to the Poisson distribution. 

It is easy to calculate the first few partial correlation 
functions by direct integration of Eq.(73): for instance, 
the probability c(2, t) is given by 


cM = (78) 


To obtain a general expression of c(n,t), we first take 
the Laplace transforms of Eq.(73) 

poo 

c(n,u ) = / dtc(l,t)e~ ut m(n,t) (79) 

Jo 

where ?n(n, t) is auxiliary function given by 

pOO 

m(n,t) = / dt'Ae- ((u+A)t '- Amax(T - t ’ 0)) TO(n-l,t') 

J max(r—1,0) 

(80) 

The initial condition is obviously, m(l,t) = 1. 

Let us introduce the generating function G m (z, t, u ) of 
the auxiliary functions m(n,t) 


G m (z,t,u) = ^ z n 1 m(n, t) (81) 

n= 1 


Inserting Eq.(80) in Eq.(81), we obtain 


pOO 

G m (z,t,u) = 1 + z / dt'G m (z,t',u) 

J max(r—£,0) 

— ((u-\-X)t' — A max(r—£,0)) 


(82) 


For t > r the generating function is constant, 
G m (z,t,y) = G m (z,T,u). For t < t, by taking the two 
partial derivatives of the integral equation Eq.(82), one 
obtains 
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d2Gm Q 2 f ' U) = z\ue-«T-» (G m (z, T-t,u) 


dt 

Simplifying we obtain 
d 2 G m (z, t, u) dG m (z,t,u ) 


^ dG m (z,T-t,u) ^ y dG(z,t,u) 


dt 2 


dt 


dt 


+ (uX + A 2 — (Az) 2 e -UT )x 


G m {z, t, u) - uX - A 2 - X 2 ze- u(T ~ t) 


(84) 


with boundary conditions (from Eq.(82)). 


y Ae 


[ Gm (Zj 0? — 1 H - 5 " u _|_^ 

I aGm ( z ’ t ’ u ' ) = zXG m (z,0,u) - X[G m (z,T,u) - 1] 

' t—T 

(85) 

whose solution is given by 


G m (z, t, u) =Ai(z, u)e Slt + Bi(z, u)e S2t 
(uX + X 2 + X 2 ze^ T ~ t '>) 


( 86 ) 


■uA + A 2 — (Az) 2 e _ 
where Si >2 are the roots of the characteristic equation 

si ,2 = ^(u ± \/(u + 2A ) 2 — 4z 2 A 2 e _ “ r ) (87) 

Finally we have 

G m (z,t,u) = (Ai(z,w)e Slt + .Bi(z,u)e S2t 
{uX + X 2 + X 2 ze~ u ^- t '>) 


uX + A 2 — (Az) 2 e _ 
+ G m {z,T,u)d(t - T ) 


0(t - *) 


( 88 ) 


where A(z,w) and B(z,u) are determined by the bound¬ 
ary conditions, Eq.(85). 

Using G m (z,t,u ) and Eq.(60) we obtain the Laplace 
transform of the correlation function. 


C(u) = / dtc(l,t)G TO (l,t,u)e ut 

Jo 

= l + \ T dtG m (l,t,u)e~ ut + G m (l,r,u)- 


u + X 
(89) 


By inserting Eq.( 88 ) in Eq.(89) we obtain 


<?(«) = 


A 


(it + A)(l + At) [ 
+ -Bi(Uw) f— 


^1(1,^) ( - 


e S2T (A + si) — u — X 


S2 


r (A + s 2 ) - u — X 


si 


+ 


(u + A ) 2 + e ut X{u 2 t + X(ut — 1)) 
u(u + A — Xe ~ UT ) 


(90) 


Figure 9(b) displays the correlation function C(t) for 
TV = 3 versus time (with t = 1). As expected for 
t < t, C(t) is constant and is equal to 1 _ | a At , because 
c(n, t) = 0, n > 1, and c(l, t) is given by Eq.(76), which is 
constant and different from 0 in this time interval. One 
also observes a discontinuity at t = t and a long time 
limit equal to 2 + 2 At+a^t 2 • We verify that, as for TV = 2, 
this is equal to 1/t with t given by Eq.(77). 

Comparing the correlation functions for TV = 2 and 
TV = 3 for the same values of A we note that the 
steady state values are higher for TV = 3 corresponding 
to a shorter time interval between particles in the 
steady state. The oscillations are more pronounced 
for TV = 2 due to the greater constraint imposed by 
the channel for smaller TV and hence greater correlations. 


VII. DISCUSSION 


The results presented in this article generalize the 
blocking model studied by Gabrielli et al. [18, 19]. In 
order to examine the situation in which blockage is trig¬ 
gered by the simultaneous presence of TV > 2 particles 
in the channel and where the particle ingress follows a 
general distribution of entry times, we have introduced 
an integral representation of the n particle survival prob¬ 
abilities. For TV = 3, we have presented exact solutions 
for the mean time to blockage, Eqs. (44,45), as well as 
the correlation functions, fluxes and other functions, for 
particles entering according to a Poisson distribution. 
For TV > 4 obtaining an exact solution appears to be 
very challenging, but we have analyzed the generic fea¬ 
tures of the model using numerical simulation. We also 
showed analytically that the mean time to blockage for 
small intensity and arbitrary TV diverges as a power of TV, 
Eq.(58). The is the result of the fact that as TV increases, 
the channel exerts a weaker constraint on the incoming 
stream and blocking is less likely. 

Future directions include the development of a mul¬ 
tichannel model, which can be applicable to filtration 
phenomenon [ 1 ], and to consider systems with diffusive 
motion that are relevant for transport through biological 
, or synthetic nanotubes [31]. 
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